program main
  use precision,only:p_
  use constants,only: m,pi   ! m is the number of momentum grid number.
  implicit none
  integer:: i,j,nstep,k,flag
  real(p_):: dp,dt,zi,theta,vte_c ! vte_c is the ratio of thermal velocity to light velocity.
  real(p_):: p(m),gamma(m),v(m),fm(m),q0(m),q(m)
  real(p_):: aa(m),bb(m),cc(m),col(m),dpp(m),djj(m)
  real(p_):: a(m),b(m),c(m),r(m)
  real(p_):: bessk,sum,k2exp
  real(p_):: w_pte,w_p0,dt_in_thermal,dt_in_p0
  real(p_):: qp(m),jp_cyclo(m),jp_landau(m)
  real(p_):: g(m),gp(m),jp3(m),sum1,sum2,pmin,vp,c10(m)
  real(p_):: q_old(m),det_q(m),det
  real(p_),parameter:: epsilon=1.0e-12
  integer,parameter:: mm=200
  external bessk
  ! this namelist provides the parameter Zi and the step number in time.
  namelist /input/ nstep,zi,theta,k2exp
  open(40,file='input_file')
  read(40,input)
  !create grid and give relativistic maxwellian distribution 
  w_pte=15._p_    !width of momentum space in unit of pte
  w_p0=w_pte*sqrt(theta) !width of momentum space in unit of p0=mc
  !w_p0=5.0_p_
  !if(w_p0<7.018) w_p0=7.018
  !write(*,*) w_p0
  dp=w_p0/(m-1)
  do j=1,m
     p(j)=0.0_p_+(j-1)*dp
     gamma(j)=sqrt(1._p_+p(j)**2)
     v(j)=p(j)/gamma(j)
     !fm(j)=1./(4.*pi*theta*bessk(2,1._p_/theta))*exp(-gamma(j)/theta)
     !fm(j)=1./(4.*pi*theta*bessk(2,1._p_/theta)*exp(1._p_/theta))*exp((1._p_-gamma(j))/theta)
     fm(j)=1./(4.*pi*theta*k2exp)*exp((1._p_-gamma(j))/theta)
     q0(j)=1.0  !initial value for q(j)
  enddo
  !calculate the difusion coefficient for maxwellian background.
  call get_diffusion_coefficient(dp,p,v,fm,dpp,djj)
  ! calculate the coefficients of the differential terms in adjoint equation
  call get_coefficient(theta,zi,dp,p,v,fm,dpp,djj,aa,bb,cc)

  !time step size in thermal collision time ( that is p0=pte )
  dt_in_thermal=1000._p_
  ! corresponding time in p0=mc collision time
  dt_in_p0=dt_in_thermal*(sqrt(theta))**3
  dt=dt_in_p0
  q=q0 ! initial condition for q

  do i=1,nstep ! advance in time
     ! calculate the collision term between maxwellian and first Legendre background.
     call get_col(theta,dp,p,v,gamma,fm,q,col)
     ! constuct the tridiagonal matrix
     call get_matrix(dt,dp,v,q,aa,bb,cc,col,a,b,c,r)
     q_old=q  !record the old value of q
     ! sovle the tridiagonal matrix equation
     call tridag(a,b,c,r,q,m)
     det_q=q-q_old  !the difference between the new and old values of q
     det=0._p_
     do j=1,m
        det=det+abs(det_q(j))  !total absolute difference.
     enddo
     if(det.le.epsilon)  then  !check whether the specified precision has been arrived at.
        flag=0
        exit
     endif
  enddo
 if(flag==1)   write(*,*) '****Not arrive at the desired precision*****'
 write(*,*) 'iteration number=',i, 'det=',det

  !Output the Spitzer function
  open(30,file='spitzer.txt')  
  do j=1,m
     write(30,*) p(j),q(j)
  enddo
  close(30)

  !calculate the derivative of q(j)
  do j=1,m-1
     qp(j)=(q(j+1)-q(j))/dp
  enddo

  do j=2,m-1
     g(j)=q(j)/p(j)   ! calculate G, defined by q/p
     gp(j)= qp(j)/p(j)-q(j)/p(j)**2   ! G prime, gp=dG/dp
  enddo

  !calculate the current drive efficiency for Landau-damped waves Eq.(18a) in Karney's paper
  ![POF 28,116, 1985]
  do j=2,m-1
     !jp_landau(j)=qp(j)/v(j)
     jp_landau(j)=(p(j)*gp(j)+g(j))/v(j)
  enddo
  !calculate the current drive efficiency for cyclotron-damped waves Eq.(18b) in Karney's paper.
  do j=2,m-1
     !jp_cyclo(j)=p(j)/v(j)*(p(j)*qp(j)-q(j))/p(j)**2
     jp_cyclo(j)=p(j)*gp(j)/v(j)
  enddo
  open(40,file='eff.txt')
  do j=2,m-1
     write(40,*) p(j), jp_landau(j),jp_cyclo(j)
  enddo
  close(40)

  ! For Eq.(19) in Karney's paper[POF 28,116, 1985]
  open(41,file='eff2.txt')
  do i=1,mm
     vp=0.01+0.98/(mm-1)*(i-1)
     pmin=vp/sqrt(1.-vp**2)
     k=pmin/dp+1
     sum1=0.0
     sum2=0.0
     do j=k,m-1
        sum1=sum1+(gp(j)*gamma(j)**2*vp**2/p(j)+g(j))*gamma(j)*fm(j)*p(j)*dp
        sum2=sum2+gamma(j)*fm(j)*p(j)*dp
     enddo
     jp3(i)=sum1/sum2/vp
     write(41,*)  vp,jp3(i)
  enddo
  close(41)
end program main


subroutine get_diffusion_coefficient(dp,p,v,fb,dpp,djj)
  ! This routine calculates the difusion coefficient for collision off isotropic background fb.
  use precision,only:p_
  use constants,only: m,pi
  implicit none
  real(p_),intent(in):: dp,p(m),v(m),fb(m)
  real(p_),intent(out)::dpp(m),djj(m)
  integer::i,j
  real(p_):: sum1,sum2,tmp

  do i=2,m-1
     sum1=0.
     do j=2,i
        sum1=sum1+p(j)**2*fb(j)*v(j)**2/v(i)**3*dp
     enddo
     sum2=0.
     do j=i+1,m
        sum2=sum2+p(j)**2*fb(j)/v(j)*dp
     enddo
     dpp(i)=4.*pi/3.*(sum1+sum2)
  enddo

  do i=2,m-1
     sum1=0.
     do j=2,i
        tmp=(3*v(i)**2-v(j)**2)/(2*v(i)**3)
        sum1=sum1+p(j)**2*fb(j)*tmp*dp
     enddo
     sum2=0.
     do j=i+1,m
        sum2=sum2+p(j)**2*fb(j)/v(j)*dp
     enddo
     djj(i)=4.*pi/3*(sum1+sum2)
  enddo

end subroutine get_diffusion_coefficient

subroutine get_coefficient(theta,zi,dp,p,v,fm,dpp,djj,aa,bb,cc)
  !This routine calculate the coefficients of the differential terms in adjoint equation
  !aa(m) for coefficient before second derivative
  !bb(m) for coefficient before first derivative
  !cc(m) for coefficient before zero derivative
  use precision,only:p_
  use constants,only: m,pi
  implicit none
  real(p_),intent(in):: theta,zi,dp,p(m),v(m),fm(m),dpp(m),djj(m)
  real(p_),intent(out)::aa(m),bb(m),cc(m)
  integer::i,j
  real(p_):: sum1,Dppdp,dvdp

  do i=2,m-1
     aa(i)=dpp(i)
  enddo

  do i=2,m-1
     sum1=0.
     do j=2,i
        sum1=sum1+p(j)**2*fm(j)*v(j)**2*dp
     enddo
     dvdp=1./sqrt((1+p(i)**2)**3)
     Dppdp=4*pi/3.*(-3./v(i)**4*dvdp*sum1)
     bb(i)=2.*dpp(i)/p(i)+Dppdp-v(i)/theta*dpp(i)
  enddo

  do i=2,m-1
      cc(i)=-2.*djj(i)/p(i)**2-zi/(v(i)*p(i)**2)
  enddo

end subroutine get_coefficient


subroutine get_col(theta,dp,p,v,gamma,fm,q,col)
  !This routine calculates the collision term between 
  !maxwellian and frst Legendre harmonics background
  !This is an implement of Eq.(7) in Karney's paper[Phys. Fluids,28(1),116,1985]
  use precision,only:p_
  use constants,only: m,pi
  implicit none
  real(p_), intent(in):: theta,dp,p(m),v(m),gamma(m),fm(m),q(m)
  real(p_),intent(out):: col(m)
  real(p_):: brace, factor1,factor2,exp1,exp2,sum1,sum2
  integer:: i,j
  do i=2,m-1
     sum1=0.
     do j=2,i
        factor1=gamma(i)/p(i)**2*   v(j)/gamma(j)**3
        factor2=gamma(i)**2/p(i)**2*v(j)/gamma(j)**3
        exp1=theta*(4*gamma(j)+6.)-1./3*(4*gamma(j)**3-9.*gamma(j))
        exp2=v(j)**2/theta*gamma(j)**3-1./3*(4.*gamma(j)**2+6.)
        brace=factor1*exp1+factor2*exp2
        sum1=sum1+p(j)**2*fm(j)*q(j)/theta*brace*dp
     enddo
     sum2=0.
     do j=i+1,m
        factor1=gamma(j)/p(j)**2*   v(i)/gamma(i)**3
        factor2=gamma(j)**2/p(j)**2*v(i)/gamma(i)**3
        exp1=theta*(4*gamma(i)+6.)-1./3*(4.*gamma(i)**3-9.*gamma(i))
        exp2=v(i)**2/theta*gamma(i)**3-1./3*(4.*gamma(i)**2+6.)
        brace=factor1*exp1+factor2*exp2
        sum2=sum2+p(j)**2*fm(j)*q(j)/theta*brace*dp
     enddo
     col(i)=4*pi*(fm(i)*q(i)/gamma(i)+0.2*(sum1+sum2))
  enddo
end subroutine get_col


subroutine get_matrix(dt,dp,v,q,aa,bb,cc,col,a,b,c,r)
  !This routine is to construct the tridiagonal matrix.
  use precision,only:p_
  use constants,only: m
  implicit none
  real(p_),intent(in):: dt,dp,v(m),q(m),aa(m),bb(m),cc(m),col(m)
  real(p_),intent(out):: a(m),b(m),c(m),r(m)
  integer:: i

  ! the following loop is to calculate 3 diagonal lines in the tridiagonal matrix
  do i=2,m-1
     a(i)=-aa(i)*dt/dp**2+bb(i)*dt/(2.*dp)
     b(i)=1.+2*aa(i)*dt/dp**2-cc(i)*dt
     c(i)=-aa(i)*dt/dp**2-bb(i)*dt/(2.*dp)
     r(i)=q(i)+dt*(col(i)+v(i))    !r(i)=q(i)+dt*col(i)-dt*v(i)/theta,for electrical conductivity
  enddo
  ! boundary conditions:
  ! the following 3 lines are to implement q(1)=0
  b(1)=1._p_
  c(1)=0._p_
  r(1)=0._p_
  ! the following 3 lines are to implement q''(m-1)=0
  a(m)=-b(m-1)-2*a(m-1)
  b(m)=a(m-1)-c(m-1)
  r(m)=-r(m-1)
end subroutine get_matrix


SUBROUTINE tridag(a,b,c,r,u,n)
  !This routine is to solve the tridiagonal matrix. This routine is from Numerical Recipes, 
  !with some minor modification.
  !(C) Copr. 1986-92 Numerical Recipes Software ,4-#.
  use precision,only:p_
  implicit none
  !  integer,PARAMETER:: NMAX=5000
  INTEGER,intent(in):: n
  REAL(p_),intent(in):: a(n),b(n),c(n),r(n)
  real(p_),intent(out)::u(n)
  INTEGER j
  ! REAL(p_) bet,gam(NMAX)
  REAL(p_) bet,gam(n)
  if(b(1).eq.0.) stop 'tridag: The tridiangle matrix equation need rewriting'
  bet=b(1)
  u(1)=r(1)/bet
  do  j=2,n
     gam(j)=c(j-1)/bet
     bet=b(j)-a(j)*gam(j)
     if(bet.eq.0.) stop 'tridag failed'
     u(j)=(r(j)-a(j)*u(j-1))/bet
  enddo
  do  j=n-1,1,-1
     u(j)=u(j)-gam(j+1)*u(j+1)
  enddo
END SUBROUTINE tridag


subroutine get_c10(theta,dp,p,v,gamma,fm,q,c10)
  !This routine calculates the collision term between 
  !the frst Legendre harmonics and Maxwellian background, c10=C(f1,fm)/(fm*cos(thet))*t0
  ! where f1=q*fm*cos(theta), t0=1/nu0, nu0 is the thermal collision frequency for electrons.
  use precision,only:p_
  use constants,only: m,pi
  implicit none
  real(p_),intent(in):: theta,dp,p(m),v(m),gamma(m),fm(m),q(m)
  real(p_),intent(out)::c10(m)
  integer::i,j
  real(p_):: sum1,qp(m),qpp(m),bb(m),dpp(m),djj(m)
  real(p_):: dvdp,dppdp
  call get_diffusion_coefficient(dp,p,v,fm,dpp,djj) !Get the diffusion coefficient for Maxwellian backgound.
  do i=2,m-1
     sum1=0.
     do j=2,i
        sum1=sum1+p(j)**2*fm(j)*v(j)**2*dp
     enddo
     dvdp=1./sqrt((1+p(i)**2)**3)
     Dppdp=4*pi/3.*(-3./v(i)**4*dvdp*sum1)
     bb(i)=2.*dpp(i)/p(i)+Dppdp-v(i)/theta*dpp(i)
  enddo

  do i=2,m-1
     qp(i)=(q(i+1)-q(i))/dp  !First derivative
     qpp(i)=(q(i-1)-2*q(i)+q(i+1))/dp**2 !Second derivative
     c10(i)=dpp(i)*qpp(i)+bb(i)*qp(i)-2.*djj(i)/p(i)**2*q(i)
  enddo
end subroutine get_c10



